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Abstract 

We develop a free energy functional for an inhomogeneous system that contains both symmetry 
conserved and symmetry broken parts of the direct pair correlation function. These correlation 
functions are found by solving the Ornstein- Zernike equation with the Percus-Yevick closure 
relation. The method developed here gives the pair correlation functions in the ordered phase 
with features that agree well with the results found by computer simulations. The theory predicts 
accurately the isotropic-nematic transition in a system of anisotropic molecules and can easily be 
extended to study other ordered phases such as smectics and crystalline solids. 
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The freezing of a fluid of anisotropic molecules into a uniaxial nematic phase is a typical 
example of a first-order phase transition in which the continuous symmetry of the isotropic 
phase is brokenfl]. In a nematic phase, molecules are aligned along a particular but arbi- 
trary direction so as to have a long range order in orientation while the translational degrees 
of freedom remain disordered as in the isotropic fluid. At the isotropic-nematic transition 
the isotropy of the space is spontaneously broken and, as a consequence, the correlations 
in the distribution of molecules loose their rotational invariance. By computer simula- 
tions of a system of ellipsoids Phoung and Schmid [2] have recently evaluated the effect of 
breaking of rotational symmetry on pair correlation functions (PCFs) and showed that in 
the nematic phase there are two qualitatively different contributions : one that preserves 
rotational invariance and other that breaks it and vanishes in the isotropic phase. The sym- 
metry preserving part of the PCFs passes smoothly without any abrupt change through the 
isotropic-nematic transition. This we expect to be true in most of the symmetry breaking 
first-order transitions. 

The correlation functions that describe the distribution of molecules in a classical fluid 
can be given as the simultaneous solution of an integral equation, the Ornstein Zernike 
(OZ) equation, and an algebraic closure relation that relates correlation functions to the 
pair potential. Well known approximations to the closure relation are the Percus-Yevick 
(PY) relation, the hypernetted-chain (HNC) relation and the mean spherical approximation 
(MSA) [3]. This approach has been used quite successfully to describe the structure of 
isotropic fluids. However, the application of the theory to ordered phases which can be 
regarded as an inhomogeneous "fluid", has so for been very limited though no features of 
the theory inherently prevent it from being used to describe structures of the ordered phases 
like crystalline solids, liquid crystals, etc. Holovko and Sokolovska [4] have used the MSA 
and Lovett equation [5] which relates one-particle density to the PCFs to solve analytically 
the OZ equation for a model of spherical particles with the long range anisotropic interaction 
and determined the PCFs in the nematic phase. However, when Phoung and Schmid[2] used 
the PY closure relation and the Lovett equation and solved the OZ equation numerically 
for system of soft ellipsoids nematic phase was not found and for this the PY closure was 
blamed. Here we adopt a method based on density-functional formalism and show that 
though the PY relation may not be very accurate it, however, leads to nematic phase with 
PCFs harmonic coefficients having features similar to those found by simulations [2] and by 
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analytical solution[4]. 

A density-functional theory (DFT) of freezing transition requires an expression of the 
grand thermodynamic potential of the system in terms of one- and two particle distribution 
functions and a relation that relates the one-particle density distribution p(x) to the PCFs. 
Such a relation is found by minimizing the grand thermodynamic potential with respect 
to p(x) with appropriate constraints [6]. The correlation functions that appear in these 
equations are of the ordered phase and are functional of p(x). The free energy functional 
that exist in the literature [7, 8] and have been used to study the freezing transitions and 
other properties of the ordered phases replace these correlations by that of the isotropic 
fluids (see ref.[6] for details). This approximation limits the applicability of the theory and 
needs improvement. 

Though in this letter we consider the nematic phase, the method developed here is equally 
applicable to other phases such as smectics and crystalline solids. We first show how the 
fact that the PCFs in the ordered phase have two distinct parts, leads us to divide the OZ 
equation and the closure relation into two sets of equations. The solutions of these equations 
give us both the symmetry conserving and symmetry breaking parts of the PCFs. Using 
these correlation functions we construct a free energy functional. This free energy func- 
tional is then used to calculate the transition parameters of the isotropic-nematic transition 
in a model system of elongated rigid molecules interacting via the Gay-Berne (GB) pair 
potential[9]. The GB potential between a pair of molecules (i, j) is written as 

«(f y ,e,,ej) = 4e(r ij ,e i ,e j )(iT 12 -ir 6 ) (1) 

where R = (r^ — <r(f y, ej) + o"o)/oo and e« is the unit vector specifying the axis of 
symmetry of the ith molecule. The expressions for the angle dependent range parameter a 
and potential well depth function e contain four parameters xq, k' , p and v. These parameters 
measure the anisotropy in the repulsive and attractive forces. The distance and energy are 
scaled by parameters a and e , respectively. The values of the parameters x ,k',fi and v 
are taken 3.0, 5.0, 2.0 and 1.0 respectively. This model system is found to exhibit first-order 
isotropic-nematic transition [ 10] . 

The OZ equation which relates the total PCF, /i(xi,x 2 ), with DPCF, c(xi,x 2 ), in an 
inhomogeneous system is written as [6,11] 

/i(xi,x 2 ) = c(xi,x 2 ) + J c(xi,x 3 )p(x 3 )/i(x3,x 2 )(ix3 (2) 
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where Xi indicates both position ri and orientation fii of the ith molecule, dx 3 = dr 3 dQ, 3 
and p(x 3 ) is the single particle density distribution. The PY relation is expressed as 

c(xi,x 2 ) = (e(xi,x 2 ) - 1)[1 + /i(x!,x 2 ) - c(xi,x 2 )] (3) 

where e(x!,x 2 ) = exp[— / 9u(x 1 , x 2 )]. The pair correlation functions h and c are assumed 
to be sum of two parts; one which corresponds to rotationally invariant and other which 
corresponds to the breaking of the rotational symmetry. Thus 

h = h {0) + h {n) and c = c (0) +c (ri) (4) 

This allows us to write the OZ and PY equations as 

hW( Xl ,X a ) = c( )(x 1 ,X 2 )+p |c(°)(x 1 ,X3) 

x/i(°)(x 3 ,x 2 )dx 3 (5) 
c( )( Xl ,x 2 ) = (e(x 1 ,x 2 )-l)[l + /i(°)(x 1 ,x 2 )- 

c (0) (xi,x 2 )] (6) 

and 

/iW( Xl ,x 2 ) = cW( Xl) x 2 ) 

+ J C (0) (x 1 ,X3)p n (x 3 )/l (0) (x3, x 2 )dx 3 

+ |(p + Pn(x 3 ))[c (0) (x 1 ,X 3 )^(x 3 ,X 2 ) 

+ C ^(x 1 ,x 3 )/i(°)(x 3 ,x 2 ) 

+cW(x 1 ,x 3 )/ i W(x 3 ,x 2 )] ( fx3 (7) 
C W( Xl ,x 2 ) = (e(x 1 ,x 2 )-l)[/ i W(x 1 ,x 2 )- 

c (ri) (xi,x 2 )] (8) 

Here p(x 3 ) = (p + p„(x 3 )) where p is the bulk number density and p„(x 3 ) = p (/(fi 3 ) — 
1). /(O) is the single particle orientation distribution function normalized to unity, i.e. 

//(n)d(n) = i. 

Eqs(5) and (6) give relations that are identical to the one used in calculating the PCFs 
in an isotropic phase. Relations given by (7) and (8) are new and as shown below give 
PCFs arising due to breaking of symmetry. To solve these equations we chose a coordinate 
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frame where the z-axis points in the direction of director h (director frame). All orientation 
dependent functions are expanded in spherical harmonics Yj m (f2)[2,12]. This yields (for 
uniaxial nematic phase of axially symmetric molecules) 

m) = 4= E /«n) (9) 

v 4tt /(cvcn) 

and 

ll2lmim,2in 

W^)^(r) (10) 



where fi = y (2/ + 1)P; and ^ stands for /i, c or e. p is the order parameter; its value is zero 
in the isotropic phase and non-zero in the nematic phase. In uniaxial symmetric phases, 
only real coefficients with m\ + m 2 — m = and even l ± + l 2 + l enter the expansion. Since 
the molecules in the model system under consideration have axial symmetry, every single / 
is even as well. Because h^°\ and e preserve the rotational symmetry, for them 

iphhimimmir) = ^i 1 i 2 i(r)C g (l 1 l 2 lm 1 m 2 m) (11) 

where C g is the Clebsch-Gordan coefficient. 

We have solved (5) and (6) for the model system of (1) using a method described in ref.[13] 
and determined the values of c[°j Jm m m (r) and h^J lm m m (r) for values of /, k up to l max = 8 
at reduced temperature T*(= ksT /e ) = 1.0 and for densities < p*(= pcr 3 ) < 0.36. 

To solve (7)-(8) we first set up linear equations for /i, ( [J a , roimaro (r) and c^ mim2m (r) using 
the expansions of eqn.(9) and (10). In these equations h\% mim2Jn (r) , c[ 1 , ) 2 , roimaro (r) and the 
order parameters P/ appear. Here we restrict ourselves to only one order parameter P 2 and 
solve these equations for < P 2 < 0.70 at the interval of AP 2 = 0.05 for all densities between 
and 0.34 and for up to l max = 4. Though we followed the same iterative method as used 
in the case of isotropic, but examined the behavior of each harmonic coefficients lm m ^ m 
at large distance and ensured that the proper convergence takes place (the computational 
details will be given elsewhere [14]). The PCFs thus generated are used in constructing the 
free energy functional. 

The reduced free energy A[p] of an inhomogeneous system is a functional of density p(x) 
and is written as [6] 

A[p] = A td [p] + A ex [p] (12) 



The ideal gas part A id [p] is exactly known, 

A ld [p] = J dxp(x)[ln{p(x)A} - 1] (13) 

where A is the cube of the thermal wavelength associated with a molecule. The excess part 
arising due to intermolecular interaction is related with the DPCF of the system as 

X 1 A 

G-X / r ~\ \ 

— — — - — r = -c xi,x 2 ; p 
dp(xi)5p(x 2 ) 

- C W( Xl ,x 2 ;[p]) (14) 

A ex [p] is found by functional integration of (14). In this integration the system is taken 
from some initial density to the final density p(x) along a path in the density space, the 
result is independent of the path of the integration[15]. For the symmetry conserving part 
c(°) the integration in density space is done taking isotropic fluid of density p«(the density 
of coexistence fluid) as reference. This leads to 

4°J[p] = A ex (pi) rfx 2 Ap( Xl )Ap(x 2 ) 

xc(xx,x 2 ) (15) 

where 

c(xi,x 2 ) = 2 J d\\ j dA'c (0) {xi,x 2 ;p2 + AA'(p - pi)} 

Ap(x) = p(x) — p h A ex (pi) is the excess reduced free energy of the isotropic fluid of density 
pi and po is the average density of the ordered phase. 

In order to integrate over [p] , we characterize the density space by two parameters A 
and £ which vary from to 1. The parameter A raises density from to p as it varies from 
to 1 whereas parameter £ raises the order parameter from to P2 as it varies from to 1. 
This integration gives 

A ex[p\ = ~\j rfxi / dx 2 p(xi)p(x 2 )5(xi,x 2 ) (16) 

where 

c(x!,x 2 ) = 4 C d& C dg C dXX f 1 d\' 
Jo Jo Jo Jo 

xc(")(x 1 ,x 2 ,AA'p ;a'P 2 ) 



It is important to note that while integrating over A the order parameter P 2 is kept fixed 
and while integrating over £ the density is kept fixed. The result does not depend whether 
integration is done first over A or £. The free energy functional of an ordered phase is the 
sum of (13), (15) and (16). Note that the Ramakrishnan and Yousuff (RY) [7] free energy 
functional is the sum of only (13) and (15) and contains an additional approximation in 
which c(xx,x 2 ) in (15) is replaced by c(x 1? x 2 ;pi). 

The grand thermodynamic potential defined as — W = A — j3p J dxp(x) where p is the 
chemical potential, is preferred to locate the freezing transition as it ensures the pressure 
and chemical potential of the two phases remain equal at the transition. The transition 
point is determined by the condition AW = W t — W = 0. The order parameters are 
determined from equations found by minimizing the grand thermodynamic potential with 
appropriate constraints [6, 14]. The isotropic-nematic transition at T* = 1.0 with one 
order parameter is found to take place at p*(— Pi^l) = 0.3325 with change in density 
Ap* (= (p — Pi) I Pi) = 0.0086 and order parameter P 2 = 0.40. If one uses the RY free 
energy functional then the transition takes place at p* = 0.3570 with Ap* = 0.0055 and 
P 2 = 0.439. The symmetry breaking part of PCFs makes the isotropic phase unstable and 
induces the emergence of the ordered phase at lower density. 

In Figs 1 and 2 we show some harmonic coefficients of DPCF in the director space 
for T* = 1.0, p* = 0.3361 and P 2 = 0.44. While the harmonic coefficient c 220 ooo(^*) 
and C22oi-io( r *) shown in Fig 1 survive both in the isotropic (P 2 = 0)and the nematic 
phases(P 2 7^ 0), the harmonic coefficients c 2 ooooo( r *) an d c o2ooo( r *) survive only in the ne- 
matic phase and vanish in the isotropic phase. The contribution arising due to symmetry 
breaking to the harmonic coefficients C22oooo( r *) an d c 2 2oi-io( r *) ar e shown by dot-dashed 
line in Fig 1 and are found to be small compared to the symmetry conserving part. Few 
selected harmonic coefficients of h are shown in Figs 3 and 4 in director space. In Fig 3 we 
plot the harmonic coefficients h 2Q0Q Q Q (r*) and /ioo2ooo( r *) which survive only in the nematic 
phase and note the oscillatory behavior which continues to survive for large values of the 
intermolecular separation r*(= r/a ). In Fig 4 we plot harmonic coefficient /i 2 2oi-io( r *) 
which is of fundamental importance and defines nematic elastic constants[l]. As shown in 
the figure it decays as l/r* at large distance. This long-range tail behaviour is attributed 
to the director transverse fluctuations [4]. This behaviour has been shown analytically by 
Holovko and Sokolovska[4] and by computer simulation by Phoung and Schmid[2]. 
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The density functional approach allows one to include more order parameters in the 
theory even though they are not included in calculating the PCFs. This is done through 
the parametrization of p(x)(see [6]). When we take two order parameters P 2 and P 4 and 
use the free energy functional developed here the transition is found to take place at T* = 1 
with p\ = 0.317, Ap* = 0.026, P 2 = 0.644 and P 4 = 0.332. These values compare well with 
the computer simulation values, p* = 0.32, P 2 = 0.66 and P 4 = 0.29[10]. This comparison 
suggests that the effect of order parameter P 4 on the PCFs is small. However, this does 
not mean that those harmonics that do not appear in the free energy functional are not 
sensitive to the values of P 4 . Inclusion of higher order parameters in calculation of PCFs is 
straightforward though computationally demanding. 

In conclusion; we developed a method of solving the OZ equation with a closure relation 
to get both the symmetry conserving and symmetry breaking parts of PCFs. Using these 
correlation functions we constructed a free energy functional to study the freezing transition 
and other properties of the ordered phase. Since the symmetry breaking parts of PCFs have 
features of the ordered phase including its geometrical packing, the free energy functional 
proposed here will allow us to study various phenomena of the ordered phases including their 
relative stability [14]. Lastly we would like to emphasize that the theory developed here can 
easily be extended to study other ordered phases such as smectics and crystalline solids. 

We thank J. Ram for his help in computation. This work was supported by a research 
grant from DST of Govt, of India, New Delhi. 
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FIG. 1: Harmonic coefficients of the DPCF in the director frame at p* = 0.3361, T* = 1.0 and 
P2 = 0.44. The dashed line shows the contribution arising due to symmetry conserving part, 
dot-dashed line the symmetry broken part and the full line the resultant. 
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FIG. 2: Harmonic coefficients of the DPCF in the director frame at p* = 0.3361, T* = 1.0 and 
P2 = 0.44. The contributions to these coefficients arise due to the symmetry breaking part only. 
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FIG. 3: Symmetry breaking h-harmonic coefficients in the director frame. Details are same as in 
Fig 2. 
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FIG. 4: harmonic coefficient /i220i-io m the director frame at p* = 0.3361, T* = 1.0 and P2 = 0.44. 
Details are same as in Fig 1. Inset shows the plot of ^2201-10 harmonics with respect to 1/r*, the 
dotted line shows the extrapolated part. 
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